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1.0 INTRODUCTION 

BOPACE is the acronym for the Boeing P/lastic Analysis Capability for 
Engines. BOPACE was developed by Boeing/Huntsville to meet the evident 
need for an advanced thermal -elastic-plastic-creep structural analyzer. 
Although BOPACE development has been strongly influenced by the 
requirements for structural analysis of engines, in particular the 
space shuttle main engine, its capabilities have been kept quite general 
and it is applicable to many types of nonlinear structures. 

The philosophy for program development was based on the following 
requi remen ts. 

1) Analysis of very high temperature and large plastic-creep effects. 

2) Treatment of cyclic thermal and mechanical loads. 

3) Improved material constitutive theory which closely follows 
actual behavior under variable temperature conditions. 

4) A stable numerical solution approach which avoids cumulative errors. 

5) Capability for handling up to 1000 degrees of freedom with moderate 
computation cost. 

Although the finite-element method was first applied to plasticity 
in the early 1960's, and several good programs for nonlinear analysis 
have since been developed, numerous improvements were indicated in order 
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1.0 (Continued) 

to satisfy the above requirements. For example, some other available 
programs assume linear plastic hardening, accumulate errors by failing 
to satisfy equilibrium at each step, or do not completely account for 
the effects of variable temperature on the elastic and plastic relations. 
The stated requirements have been effectively met by the current BOPACE 
program version. In addition, the research and development effort has 
led to an improved hardening theory for cyclic plasticity, a method for 
representing general cases of load reversal, and advanced techniques for 
improving the accuracy and controlling convergence of highly nonlinear 
sol utions. 

Two versions of the current BOPACE program are available. The first 
is a 300-D0F version developed for fast analysis of small size problems 
within moderate core-storage limitations. The second is the basic 1000- 
DOF version. In addition, a low-core modification of the 1000-DDF 
version has been accomplished through the use of overlays and dynamic 
storage of arrays. BOPACE is written in FORTRAN IV and has been 
extensively run on both the IBM 360 and IINIVAC 1108 computer systems. 
Documentation consists of three volumes: Theoretical Manual, User Manual 

(including example problems), and Programmer Manual. 

The BOPACE development and programming effort has been performed at 
Boeing/Huntsville by Dr. R. G. Vos, with suggestions and review by 
W. H. Armstrong. A. H. Spring assisted with many analyses and program 
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1.0 (Continued) 

checkout. J. L. Ballinger of Boeing Computer Services modified and 
programmed the Gauss wavefront solution method. Recognition is also 
due to N. L. Schlemmer, L. Salter and R. Hurford at the NASA Marshall 
Space Flight Center, and L. Johnston of Brown Engineering Co., for 
their suggestions and support of the program development. 
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2.0 MATERIAL CONSTITUTIVE THEORY 

The basic purpose of classical constitutive theory in an el asto-visco- 
plastic program such as BOPACE is to provide incremental relations 
between stresses and strains. BOPACE uses these relations with the finite- 
element stiffness method to provide a convenient and efficient approach 
for solution of an important class of nonlinear material problems. 

BOPACE performs plane-stress and plane-strain analyses, as well as limited 
3-dimensional analyses involving prescribed non-zero values of normal 
stress or strain. It accounts for elastic, plastic, thermal and creep 
deformations, and the nonlinear dependence of all deformations on 
temperature. New constitutive developments include a combined plastic- 
hardening theory, and a generalized approach to cyclic load reversal. The 
BOPACE constitutive theory is developed by a tensorial approach which 
provides all relations in a form which is invariant under coordinate 
transformations. 

2.1 ELASTICITY EQUATIONS 

This section defines the cumulative and incremental forms of the relations 
for temperature- dependent elasticity. The BOPACE formulations for plane 
stress and plane strain are developed, including the effects of prescribed 
non-zero values for normal stress or strain. 

General Concepts and 3-D Relations - The basic cumulative stress-strain 
relation, for either temperature-dependent or temperature-independent 
elasticity, is 
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2.1 (Continued) 


a 


ij 




(2.1-1) 


vhere o and e e are the 3x3 tensors of stress and elastic (recoverable) 
strain, respectively, and D e is the tensor of elastic coefficients 
which may depend on temperature. For convenience we will use the 
equivalent single-subscript notation 



( 2 . 1 - 2 ) 


where subscripts i and j range over all nine of the tensor components. 
For 3-dimensional analysis the relation 2.1-2 is taken as 
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2.1 (Continued) 
where 
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E 
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Here E is Young's modulus and v is Poisson's ratio. 


Note that the elasticity matrix in Equation 2.1-3 is consistent with 

the tensorial definition of shear strains (e.g. e e = y e .../2, 

xy xy 

where y is the engineering definition of shear strain). Tensorial 
xy 

definitions are used throughout the BOPACE program in order to easily 
formulate constitutive theory which is invariant with respect to 
coordinate transformations. Although such invariance can be achieved 
by a careful use of engineering strain definitions, such definitions 
are probably responsible for many invariance difficulties such as 
those described in Reference 1. 
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2.1 (Continued) 

The last three of Equations 2.1-3 are somewhat redundant and may be discarded 

to give an abbreviated 6-component form (e.g. a = a and e = e ). 

y x xy y a xy 

It should be emphasized, however, that in performing later summations 
all non-zero values of the nine tensor components must be accounted for. 


Relations for Plane Stress and Plane Strain - For plane stress, the 
abbreviated form of Equation 2.1-2 is 



and for plain strain it is 


"1 

V 

0 

/ e \ 

( G xx ) 

V 

1 

0 

) £ e l 
)yy( 

_0 

0 

l-v_ 

U e ) 

xy 


( a xx 

\ 

"1-v 

V 

0 

e 

f xx 'j 

|°yy 

{ _ E 

V 

1-v 

0 

} G f 

< e yy 

t (l+v)(l-2v) 




a *y 

) 

0 

0 

1 

ro 

1 

f e 1 

£ xy ' 


(2.1-4) 


(2.1-5) 


Incremental Relations - For the case of temperature-independent 
elasticity the incremental stress-strain relations are simply 
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2.1 (Continued) 

AOj = D®. Ac® (2.1-6) 

I l J J 

where A denotes an incremental quantity and D e is the appropriate 
elasticity matrix. 

When temperature dependence is considered, the incremental relations 
become 

pQ 

where e are the elastic strains at the beginning of the increment, 

p "J 

D is the elasticity matrix evaluated for the temperature at the 

p pi p Q 

end of the increment, and AD = D - D is the change in elasticity 
matrix from beginning to end of the increment. The first term in 
Equation 2.1-7 accounts for stress change due only to change in elastic 
properties, while the second term accounts for additional stress change 
due to the increment of elastic strain. 

BOPACE Formulations - The BOPACE program allows the option for either 
a "plane-stress" analysis with prescribed non-zero values of normal 
stress o zz , or a "plane-strain" analysis with prescribed non-zero values 
of normal strain 

In plane-strain problems, BOPACE uses the 3-D elasticity matrix of 
Equation 2.1-3 because, even though the total (physical) strain c zz 
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2.1 (Continued) 

might vanish, the elastic strain e® z is in general not zero. If a 
non-zero value is prescribed by the user for c , no additional 
difficulties are introduced. 


For plane-stress problems, BOPACE uses the elasticity matrix of 
Equation 2.1-4. For the general case when a non-zero value is prescribed 
by the user for ° zz > an additional term is required in the stress-strain 
relation. The resulting general plane-stress equation is 



E 

1-v 2 


" 1 

V 

_ 0 


V 

1 

0 



(2.1-8a) 



(2.1 -8b) 


The general incremental plane-stress relation for temperature-dependent 
elasticity then follows Equation 2.1-7: 


Acr. = AD?. e?° 

l 1 3 J 


D e ! Ac? + 
ij 3 




zz 


6 1 

+ ]u A a 


zz 


(2.1-9) 


where c zz and Ao zz are respectively the initial and incremental values 
of prescribed normal stress. 
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2.2 THERMAL STRAIN 

Alternate Formulations - The conventional description of thermal strain 
is given by 



where e 1 " denotes thermal strain, T is the temperature, and y is the thermal 
coefficient of expansion which may be a function of temperature. 


An alternate integrated description of thermal strain is 



( 2 . 2 - 2 ) 


where here e* gives the thermal strain directly as a function of temperature. 

,+ 

If only incremental thermal strains are of interest, e may be taken 
as zero at any convenient reference temperature. 


BOPACE Formulation - BOPACE uses the direct form 2.2-2. This form is 
preferred over that involving a thermal expansion coefficient because 
accumulated errors in thermal strain are not introduced. These errors 
could arise with the form 2.2-1, in case y varied with temperature and 
the specified heating and cooling sequences used different temperature 
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2.2 (Continued) 

increments. BOPACE takes the structural fabrication temperature as the 
reference temperature for zero thermal strains. 


2.3 PLASTICITY 

This section defines the incremental elasto-plasti c relations used in the 
BOPACE program. BOPACE employs a new concept of combined isotropic and 
kinematic hardening, and accounts for temperature-dependent elasto-plastic 
behavior as well as a generalized form of cyclic load reversal. In order 
to develop the constitutive theory in a straightforward manner, discussion 
of the effects of temperature-dependent elasticity on the elasto-plastic 
relations is deferred until Section 2.5. 


Definitions - The following nomenclature is defined. 



total stress 

stress center (of yield surface in kinematic 
hardening) 

deviatoric (total - hydrostatic) stress 
deviatoric stress center 
s - a = relative deviatoric stress 
elastic (recoverable) strain 

plastic (time independent non stress-inducing strain) 
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2.3 (Continued) 

General Concepts - The basic concepts in most el asto-plasti c theories are 
those of a yield surface , the dependence of yield on only the deviatori c 
stress components, incompressibility under plastic strains, and normality 
of the incremental plastic-strain vector to the yield surface. The 
definition of a particular theory requires assumptions for three basic 
constituents: 

1. a surface relating the stress components at yield 

2. a flow rule defining a direction for the incremental plastic- 
strain vector 

3. a hardening rule. 

Yield Surface - BOPACE employs the Huber-Mises yield surface [2], 
defined by the relative deviatoric stresses as 

/N /*S /\ /N 

F = s.s. - s?s° = 0 (2.3-1) 

where the summation is again taken over all nine tensor components of s. 

A 

The s? are components of a point on the yield surface at a known condition 
of temperature and plastic deformation, e.g. from a uniaxial test. 

Equation 2.3-1 holds whenever the material is plastic, i.e. whenever 
the components of s are on the yield surface. Function F may be thought 
of geometrically as defining a hypersphere in 9-dimensional stress space. 
Alternatively, when expressed in the 3-D space of principal stresses, this 
yield surface can be shown to be an open-ended circular cylinder whose 
axis passes through the origin and makes equal angles with each of the 
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2.3 (Continued) 

three principal stress axes. The Huber-Mises yield surface is generally 
used to describe plasticity in metals because it agrees reasonably well 
with test results and it gives a smooth surface which is convenient for 
calculations. 

Flow Rule - BOPACE uses the Prandtl -Reuss flow rule, which is the usual 
rule associated with the Huber-Mises yield surface. The assumptions are 
that the material is incompressible under plastic flow, and that increments 
of plastic strain are normal to the yield surface at the stress point. 

These provide the relation 

AeP = As. (2.3-2) 

where A is a flow parameter. 

Basic Hardening Concepts - An elastic-plastic material which work hardens 
in the plastic range is commonly analyzed using either of two classical 
hardening theories. Isotropic hardening [ 3 ], which assumes a uniform 
expansion of the yield surface during plastic flow, accounts for change 
in size of the hysteresis loop during cycling. Kinematic hardening [4], 
which assumes a rigid translation of the yield surface in the direction of 
the plastic strain increment, accounts for the pronounced Bauschinger effect 
which is evident in cyclic behavior of most metals. In general, the actual 
cyclic behavior can be more accurately described by a combination of 
isotropic and kinematic hardening. A combined hardening theory has. been 
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2.3 (Continued) 

given by Hodge [5] for materials which satisfy the Tresca yield condition. 
Because a better representation for most metals is provided by the Huber- 
Mises yield surface, a corresponding combined hardening theory [6] has 
been developed for the BOPACE program. 

Hardening Parameters - A simple combined hardening theory such as that 
presented in Reference 6 makes two basic assumptions: 

1) Size of the yield surface is a function of a cumulative hardening 
parameter, k. This means that the isotropic hardening, i.e. the 
incremental change in size of the yield surface, depends on the 
initial value of k and its change Ak. 

2) Yield surface transl ation is related (but only in an incremental 
manner) to a kinematic hardening parameter, k^. The kinematic 
hardening, i.e. the incremental translation of the yield surface, 
depends on the initial value of and its change Ak^. 

It will be evident in the discussion to follow that isotropic hardening 
can be related to k on either a cumulative or incremental basis, while 
kinematic hardening can be related to k ^ only on an incremental basis. 

In addition to the parameters k and k^, hardening is also a function of 
temperature. 

Figure 2.3-1 shows hysteresis loops for the first two strain-controlled 
cycles of a typical material which exhibits combined isotropic and 
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FIGURE 2.3-1 : COMBINED HARDENING BEHAVIOR (NON-STABILIZED) 


a 



FIGURE 2.3-2 : VARIABLE TEMPERATURE HYSTERESIS LOOPS (STABILIZED) 


o 



FIGURE 2.3-3 : VARIABLE TEMPERATURE HARDENING EFFECTS 


2-12 


D5-1 7266-1 


2.3 (Continued) 

kinematic hardening. The Bauschinger kinematic hardening effect is 
apparent in that the initial yielding in tension causes a reduced yield 
stress in compression, i.e. a shift of the yield center by an amount a. 
Successive yielding in compression causes a reduced yield stress in 
tension, and so forth. Isotropic hardening causes the increase in size 
of the hysteresis loop with continued cycling. The hysteresis loops for 
many materials become stabilized after a number of cycles, and then may 
begin to decrease in size as further deformation causes a softening 
effect. 

Figure 2.3-2 shows the stabilized hysteresis loops for a material at 
various temperature levels. (Different strain ranges are used to separate 
the loops for purpose of illustration). The hysteresis loop of a material 
typically decreases in size with increasing temperature. Note that the 
size of the yield surface will vary in a similar manner with temperature. 
Also the rates of isotropic and kinematic hardening with respect to 
plastic deformation vary with temperature. 

The isotropic hardening parameter k may be appropriately taken as either 
the cumulative plastic work , or as the sum of increments of effective 
plastic strain . The kinematic hardening parameter ic* 4 must account for 
the Bauschinger effect in cyclic loading, and it may be taken as ah adjusted 
value of k. As long as no load reversal occurs and the loading is 
proportional, is simply equal to k. However, must be set to zero 

at the start of each increment in which a complete load reversal occurs. 
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2.3 (Continued) 

(A complete load reversal occurs when the incremental plastic strain 
vector has a direction exactly reversed from that of the previous 
plastic increment). For an incomplete load reversal, the BOPACE program 
computes the starting value for <k by multiplying the existing accumulated 
value of by the factor (1 + COS IN E ) /2 , where COSINE is the Cosine of 
the angle between successive incremental plastic strain vectors. At the 
end of each increment, becomes + Ak. 

Because the Bauschinger effect varies with cumulative deformation in 
certain materials (e.g. it may become more pronounced as plastic cycling 
continues), BOPACE allows an additional option for the kinematic hardening 
to be defined as a product of two functions. The first is a function of 
and defines the shape of the kinematic hardening, while the second is 
an additional factor which depends on k and defines the magnitude of the 
kinematic hardening. 

In order to implement the BOPACE hardening theory, it must be determined 
how the size of the yield surface varies with temperature. In addition, 
the dependence of isotropic and kinematic hardening on the parameters k 
and must be determined. This is accomplished through cyclic tests 
performed at various temperature levels. Because the cyclic hardening 
behavior is thus determined for a constant temperature, an assumption 
must be made for variable temperature cycling. The hardening effects of 
variable temperature are illustrated in Figure 2.3-3. As long as 
temperature remains constant, plastic hardening behavior is defined by 
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2.3 (Continued) 

following the shape of a stress-strain curve at the given temperature, 
say to the point 0 on the T2 curve. If temperature changes to T3, and 
then plastic deformation continues, an initial point must be determined 
on the T 3 curve from which the new yield surface size and initial 
hardening slope may be determined. This transfer from curve T2 to 
curve T3 requires a definition of the basis for hardening, i.e. the 
definition of the parameters k and k^. BOP ACE allows the option of 
either plastic work or the sum of increments of effective plastic strain 
to be used as the hardening basis. The strain and work options correspond 
to the respective points 1 and 2 in Figure 2.3-3. 

The hardening relationship determined from a series of cyclic tests may 
depend somewhat on the strain range used in a particular test. If strain 
range is a significant factor the test conditions should duplicate the 
approximate expected strain range for which an analysis is to be made. 

The choice between plastic strain and plastic work as a basis for the 
hardening parameters k and may depend to a large extent on which basis 
provides the better overall representation of cyclic behavior at various 
strain ranges. 


Multiaxial Hardening Rule - The kinematic hardening rule employed in 
BOPACE is that due to Prager [4]. It gives the increment of yield surface 
translation in terms of the incremental plastic strains, as 


■ c u 4E j * f c 'u 4E j 


(2.3-3) 
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2.3 (Continued) 

where c is the kinematic contribution to the slope of the uniaxial 
stress vs. plastic-strain curve, and I is the identity matrix. An 
alternate hardening rule due to Ziegler [ 7 ] is preferred by some 
plasticity analysts because the form of Ziegler's rule does not change 
with reduction in the number of spatial dimensions, and it is therefore 
supposed to simplify the calculations. Prager's rule is considered 
more acceptable from a physical point of view, however, and it presents 
no difficulty when all components of the required tensors are retained 
as in the BOPACE program. Note that for Prager's kinematic hardening 
rule, the deviatoric stress center is equal to the stress center, 

a i = a i* 

The isotropic hardening, i.e. change in size of the yield surface due 
to plastic deformation, is defined for a proportional test loading by 

' R ?j ' § rI U iE j (2 - 3 - 4) 

where r is the isotropic contribution to the slope of the uniaxial 
stress vs. plastic-strain curve. 

The necessary condition that stresses remain on the yield surface is 
satisfied by taking the differential of Equation 2.3-1. The condition 
is AF = 0, which can be shown to give 

S. Ao i - S.j Aa.. - S? As? = 0 (2.3-5) 


or 
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2.3 (Continued) 


s.j Ao.j - AX = 0 


(2.3-6) 


where A = j s. Ac^ + js? As!? 


C. . s. s . + R? • s? s? 
ij i J iJ i J 


(2.3-7) 


The key to a successful combined hardening theory is the proper 
determination of the hardening variable A. The BOPACE program uses 
hardening tables which give the yield-surface size and the surface 
translation as functions of the hardening parameters k and « k . These 
are two-dimensional tables for each material whose ordinates and 
abscissas are, respectively, temperature and hardening parameter. Given 
the initial values of < and K k at the beginning of an increment, and 
estimated values for Ak and A< k , the corresponding increments of isotropic 
and kinematic stress increase are obtained from the hardening tables. 
(Hardening due to temperature change is included by adding it to the 
isotropic stress increment). The hardening slopes c and r are then 
computed, by dividing the incremental stress increases by the estimated 
increment of effective plastic strain. This procedure gives average values 
for the slopes c and r during the increment, and tends to produce a quite 
stable numerical iterative process. Note that it is the isotropic and 
kinematic stress increases, rather than the slopes c and r, which are 
directly relatable to the hardening parameters. The choice of a test value 
for s° in Equation 2.3-7 is arbitrary, as long as it is a point on a yield 

/N 

surface of size corresponding to s, i.e. a surface with equal values of 

A 

temperature and parameter k. It is convenient in BOPACE to take s' 1 equal 

A 

to s. 
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2.3 (Continued) 


-Incremental Stress-Strain Relation - The incremental stress-strain 
relation now follows the development of References 8 and 9. Take 


Ac. 


1J 


A£ 1 


= D®. 
lj 


Ae 


e+p 


- D 


1J 




(2.3-8) 


where D e is the appropriate matrix of elastic constants. 


Then 


AA = s. A a^ : 


S. D e . Ae e+ P - s. D®. S.A 
1 TJ J 1 1J J 


(2.3-9) 


which gives 

A = s i Ae® +P /(A + s k D®^ s fc ) 


(2.3-10) 


Substituting Equation 2.3-10 into Equation 2.3-8 provides the desired 
relation 


or 


A°i 


D ik s k S £ D lj 

^ P 

A + s D s n 
m mn n 


AcfP 


Ac. 


(D®. + DP.) Ae® +P = D. • Ae® 
v ir j ij j 


e+p 


(2.3-11) 


(2.3-12) 


D is the elas to-plastic Jacobian (tangent-stiffness) matrix relating 
incremental stresses to incremental elastic+plastic strains. In effect, 
it separates the elastic and plastic strains and determines the incremental 
stress corresponding to the incremental elastic strain. D p is the stiffness 
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reduction due to plastic flow, and becomes zero for the case of infinite 
hardening, i.e. A = or equivalently the total slope (c+r) of the 
stress vs. plastic-strain curve is infinite. 


Yield Surface Correction - Because an incremental plasticity solution 
generally involves a variable stress-strain slope and a changing direction 
of the surface normal within each increment, the yield surface relationship 
given by Equation 2.3-1 may not be satisfied exactly. It is probably 
desirable to exactly satisfy this relationship, i.e. have F = 0, at the 
end of each plastic increment. This can best be accomplished by making 
a correction to the stress values, and holding constant all other quantities 
involved in Equation 2.3-1. Proceeding in a manner similar to that of 
Reference 10, the change in F due to a change in stress is 

AF = 2 s. Aa i (2.3-13) 


We take the stress correction to be defined by 

Aa i = pn i 


(2.3-14) 


where p is a scalar factor to be 
while 



determined for the stress correction. 


(2. 3-1 5a) 
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2.3 (Continued) 
for plane strain, and 

xx 1 
yy ( 

) (2. 3-1 5b) 

0 \ 

xy / 

for plane stress. Thus the direction of the stress correction is taken 
normal to the yield surface, except that in plane stress the value of 

°zz is not correc ' tec1, l" hen 

AF = 2 n 1 -n i P (2.3-16) 



giving 


P 


AF 

2 n 1 n 1 


(2.3-17) 


and 


Ao^ = pn i 


AF 


2 n i n i 


(2.3-18) 


The stress correction is computed from Equation 2.3-18 using for AF the 
negative of the error evaluated for function F in Equation 2.3-1. Because 
the yield surface relation is to be satisfied at the end of the increment, 
all quantities in Equation 2.3-18 are also evaluated at the end of the 
increment. 
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2.3 (Continued) 

Effective Stress-Strain and Plastic Work - The concepts of " effecti ve 
stress " and " effective strain" are related to pi as ti c work , and are used 
extensively in the literature on plasticity. These concepts are used in 
a limited way in the development of constitutive theory for the BOPACE 
program. 

Because they can easily be misapplied, especially in the presence of 
kinematic hardening, the use and limitations of the concepts are briefly 
discussed here for the Mises plasticity theory. 

Due to characteristics of the Prager hardening theory, the following 
statements of equivalence and proportionality should first be noted. 

Aa i - AefiwSjjjfe s i (2.3-19) 

Because of the incremental nature of kinematic hardening, s^ and 
are in general not proportional. 

The Mises effective stress a is defined by 

o 2 = | s ■ s i (2.3-20) 

The incremental and cumulative values for plastic work, W p , are given by 
AW P = a. Ae p (2.3-21a) 
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2.3 (Continued) 
and 

W P = 2 a. Ae? (2.3-21b) 

where Z denotes summation over all increments. For the special case of 
proportional loading (i.e. loading in which all stress components are 
increased proportionately) followed by a constant stress level (i.e. no 
plastic hardening), the cumulative plastic work is given by 

W p = o. ej? (2. 3-21 c) 

As a matter of convenience in computing plastic work, an increment of 
effective plastic strain, Ae^, has historically been defined by 

(Ae p ) 2 = | Ae? Ae? (2.3-22) 

At this point, however, care must be exercised in using the historical 
calculation for plastic work. If kinematic hardening were zero, then 

'V _ /v 

s .j = s • , and because Ae p is proportional to s.. the use of Equations 2.3-20 
and 2.3-22 would give plastic work as 

AW P = a A? (2.3-23a) 

and 

W p = E a Ae^ (2.3-23b) 
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2.3 (Continued) 

If in addition, the condition were one of proportional loading and 
constant stress, then by defining the cumulative effective plastic strain, 
e* 3 , in the same manner as A? 5 , we would have 

W p = a? (2.3-23c) 

Of course the Equations 2.3-23 in general are not valid , because of the 
presence of kinematic hardening and non-proportional loading. Thus plastic 
work must be computed from Equation 2. 3-21 a and b, rather than from the 
product of effective stress and strain quantities. 

The quantity I 13 serves little purpose in a general plasticity analysis, 
although it is a tensorially invariant quantity and does provide a measure 
of net residual deformation. For a rational measure of deformation history, 
either the plastic work, W* 3 , or the sum of increments of effective plastic 
strain, 2 Ac* 3 , is appropriate. The difference in concept between the 
quantities W* 3 and E Ae * 3 should, however, be recognized. 

2.4 CREEP 

Stages - Metals characteristically exhibit the three stages of primary, 
secondary and tertiary creep. Figure 2.4-1 shows these stages in a typical 
creep history under conditions of constant temperature and stress. 

Because creep rate varies considerably during the different stages, the 
description of actual creep histories is considered to be essential for 
an accurate analysis. The BOPACE program accounts for the creep time 
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FIGURE 2.4-1 : TYPICAL CREEP STAGES 



FIGURE 2.4-2 : BOPACE CREEP REPRESENTATION (EXAMPLE FOR VARIABLE TEMPERATURE 

AND CONSTANT STRESS) 



FIGURE 2.4-3 : BOPACE CREEP HARDENING OPTIONS (EXAMPLE FOR CONSTANT 

TEMPERATURE AND VARIABLE STRESS) 
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2.4 (Continued) 

history by allowing the user to define, by a series of input points, 
the shape of the effective-creep-strain vs. time curve for each material. 

Temperature and Stress Effects - Creep rate in most metals is very 
dependent upon temperature and stress level . The BOPACE approach to 
creep analysis provides a reasonable description of temperature and stress 
effects, while avoiding excessive storage and computational requirements. 
For each material, BOPACE requires a reference creep curve shape which 
gives the relative variation of effective-creep-strain vs. time for the 
various stages considered. This shape is assumed to be valid for all 
the temperatures and stress levels of the particular material. A table 
of creep factors for the material is then specified as a function of 
temperature and effective stress, and a portion of the actual creep curve 
is determined by multiplying the reference creep curve by the appropriate 
factor using the average temperature and stress during the increment. 
Figure 2.4-2 shows portions of typical creep curves for the special case 
of constant stress level and variable temperature. Note that according 
to BOPACE assumptions these curves have the same shape. 

Hardening - As long as the temperature and stress level remain constant, 
an increment of creep is determined by following the corresponding creep 
curve for the given time increment. However, if temperature or stress 
level changes, an initial point must be identified on the corresponding 
new creep curve in order to determine the new creep rate. This transfer 
from one curve to another requires an assumption for creep hardening. 
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2.4 (Continued) 

which in BOPACE is defined by a single hardening parameter, k . BOPACE 

c 

allows the option of either age, strain, or work hardening, for which k 
is defined respectively as the accumulated time, sum of increments of 
effective creep strain, or creep work. These options are illustrated in 
Figure 2.4-3 for a case of constant temperature. Creep during the 
preceding increments has progressed to the point 0 on the o-| curve. The 
average effective stress during the present increment is a 2 , which gives 
the initial points 1, 2 and 3, respectively, for the options of age, strain 
and work hardening. Incremental creep for the current increment is then 
determined by continuing along the curve from the appropriate initial 
point, for a distance equal to the specified creep time increment. In 
the general case both temperature and stress will vary from one increment 
to the next, but the hardening option still determines in the same manner 
how the transfer is made between the creep curves. 

Load Reversal - The main use of the creep-hardening parameter k c 
comes into play during a load reversal. When a complete reversal occurs, 
k c is set to zero and the initial point on the creep curve is taken as that 
corresponding to a zero value of k c . (A complete load reversal occurs 
if the incremental creep-strain vector has a direction exactly reversed 
from that of the preceding creep increment.) For an incomplete load 
reversal, the BOPACE program computes the starting value for k c by 
multiplying the existing value of k c by the factor 0 + C0SINE)/2, where 
COSINE is the Cosine of the angle between successive incremental creep 
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2.4 (Continued) 

strain vectors. Parameter k c then accumulates as before, i.e. at the 
end of each increment k c becomes k c + Ak c . 

Multi axial Flow Rule - The incremental creep-strain vector has 
historically been taken normal to a Mises type of surface which passes 
through the stress point. When kinematic plastic hardening is 
considered, this surface could be taken either as the actual translated 
yield surface, or as an untranslated surface which passes through the 
stress point but whose center remains at the origin. The appropriate 
choice of surface is not clear, and because the untranslated form of 
surface is somewhat simpler, further complications are perhaps unjustified. 
BOPACE therefore defines the multi axial creep flow rule by 

AeJ = (|a^V)s. (2.4-1) 

where As 0 is the increment of effective creep strain defined by 

(Ae^ 2 = | Ae? Ac? (2.4-2) 

The increment of creep work is then given by either 

AW C = a i Ae9 (2.4-3a) 

or 

AW C = o’ Ae^ (2.4-3b) 
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2.5 COMPLETE STRESS-STRAIN RELATIONS 


In Sections 2.1 to 2.4, the basic theory used in BOPACE for elasticity, 
thermal strains, plasticity and creep has been discussed. The present 
section describes the complete stress-strain relations, and the manner 
in which simultaneous elastic, plastic, thermal and creep strains are 
accounted for. The combined effects of temperature-dependent elasticity 
and plasticity are included, and the specializations are described for 
plane-stress and plane-strain problems. 


General 3-D Relations - For temperature-dependent behavior, the required 
elasto-plastic incremental stress-strain relation follows from Equations 
2.1-7 and 2.3-8: 

■,el 


. ^ , n e eO 

Aa.j - AD^ j e . 


+ D?'. Ae® +P 
J 




(2-5-1) 


ij J 

Here the first term accounts for stress change due to change in elastic 
properties, while the second and third terms account for stress change due 
to change in elastic strain. Following Equation 2.3-9, 


Aa = s. Ac i « s. ad?, e ?° + s, D«] A E ? +P - S, D?! s, 


13 E 3 


1J 


U J 


where again 


r s s + R° s° s° 

C ij i J R ij i J 


(2.5-2) 


(2.5-3) 


For the general case of temperature-dependent plasticity, R° accounts for 
isotropic hardening due to both plastic deformation and temperature. Then 


; Anfi eO , 1 n el . e+p 

$ 1 A \l e .i * 5 i c 1,i Ae .i 

A ♦ s k B tl \ 


(2.5-4) 
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2.5 (Continued) 


Substituting Equation 2.5-4 into 2.5-1 gives 


D?J s,. s. AD® 


d:: s, s. D! 


ik a k eO + D el _ ik s k u i5 


A = AD= - - uy + D? . ~ e1 

1 1 1J A + s D e s / J \ 1J A + s D e S 

m mn n / \ m m n 


or, using abbreviated notation 

Ao. - (AD®. + A D P.) e®° + (D®j + aP])Ae® + P = 


AD. . e^° + D 1 . • Ae® +P 
1J J 1J J 


A e+ P 

AC j 


(2.5-5) 


(2.5-6) 


Thus the increment of stress is given by the sum of two products: an 
incremental matrix times the initial elastic strains, plus an end-of- 
increment matrix times the incremental elastic+plastic strains. 


Calculation of Strains - The total strain increment, Ae, is computed from 
the incremental displacements. The strain Ae is the physical (measurable) 
strain in the structure, and is the sum of the elastic, plastic, thermal 
and creep strains. The elastic+plastic strain, Ae e+P , is needed for the 
calculation of the stress increment in Equation 2.5-6. It is computed by 

Ae® +P = Ac- - Ae* - As® (2.5-7) 

t c 

where e and e are the thermal and creep strains, respectively. Because 
Ae® depends on average stress level during the increment, and because the 
normal s to the yield surface as well as the plastic slope parameters c and 
r vary somewhat during the increment, BOPACE uses an iterative approach for 
solution. The solution approach is discussed in Section 4. 
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2.5 (Continued) 


Plane Stress and Plane Strain - The development of Equations 2.5-1 through 
2.5-6 is based on 3-dimensional relations in which all nine components of 
the various tensors are properly accounted for. For plane-strain problems, 
BOPACE uses the full 3-D relations and thus accounts automatically for 
prescribed non-zero values of the strains e zz . 

For plane-stress problems, BOPACE uses the abbreviated elasticity matrix 
of Equation 2.1-4 in Equation 2.5-6 (except that both the xy and yx tensor 
components are accounted for). The hardening variable A must be computed 
from Equation 2.5-3 using all components of s. In the case of a specified 

p Q gl 

non-zero value of the normal stress a zz> the quantity Ap.. a zz + p. Aa zz 
is added to the right-hand-side of Equation 2.5-1. Following through the 
incremental stress-strain derivation as before gives 


, _ : , n e eO , : n el . e+p , ; 

A - s. AD ij e. + s. D. . ACj + Si 


APi a zz + s i u i ] Ao zz + s zz Ao zz (2 * 5_R) 


A ♦ s k D kl 


The plane-stress equivalent of Equation 2.5-6 is then obtained in a similar 
manner as before. 
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3.0 FINITE ELEMENT REPRESENTATION 

BOPACE uses the Constant-Strain-Triangle (CST) element in its thermal, 
elastic, plastic and creep analysis. The following sections describe 
the element stiffness matrices, the coordinate systems employed, and 
the various coordinate-related transformations. 


3.1 CST ELEMENT STIFFNESS MATRICES 

Element Coordinates and Displacement Functions - The CST element is 
shown in Figure 3.1-1. The element is defined using the x-y Cartesian 
coordinate system, which is arbitrarily oriented with respect to the 
element. (A coordinate system with origin at node 1, and x-axis along 
the nodal line 1-2, is used for output of element stresses and strains.) 
Displacements u and v correspond to the x and y coordinate directions, 
respectively. 


Following the general development of Reference 11, the displacement 
functions are defined by 




b 2 



(3.1-1) 


giving linear variation of displacements and constant strains within the 
element. 
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3.1 


(Continued) 



FIGURE 3.1-1: CST ELEMENT AND COORDINATES 


Strain - Displacement Matrix - Evaluating the unknown coefficients b] 
through b 6 in terms of nodal displacements, substituting these values 
into Equation 3.1-1, and differentiating to obtain expressions for the 
strains, gives 


'•xx 
£ yy 
£ xy 
£ yx 


au/ax 

9V/3y 

(3u/ay+av/3x)/2 

(au/3y+3v/3x)/2 


[r] 


U 1 

V 1 

u 2 

v 2 

u 3 

v 3 


where the strain-displacement matrix, r, is given by 

0 y 31 0 

x 32 0 " x 31 

x 32/ 2 -y32/ 2 " x 3l/ 2 y3l/ 2 x 2 l /2 

x 32/2 -y 32 /2 -x 31 /2 y 3]/2 x 21 /2 


[r] = — 
L J 2A 


-y 32 

0 


-y2i 

0 


( 2 . 1 - 2 ) 


0 

X 21 

-y2i/2 

-y 2] /2 
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3.1 (Continued) 

Here A is the area of the triangle, x-jj and y-j j are abbreviated notation 
for the differences x-j-xj and y.j-yj, respectively, and and are 
displacements at the ith node in the x and y directions, respectively. 


Plane-Stress Stiffness Matrix - For plane-stress conditions it can be 
shown that the CST stiffness matrix, k, is given by 


k ij A D mn r mi r mj 


The elasticity matrix D is 


[D] = 


1-vZ 


1 v 0 0 

v 1 0 0 

0 0 1-v 0 

00 0 1-v 


(3.1-3) 


where E is Young's modulus and v is Poisson's ratio- 


The stiffness matrix defines element generalized forces, p^ , in terms of 
displacements, q n - , by 

Pi = k i j (3.1-4) 

where q is the right-hand-side nodal -displacement vector in Equation 
3.1-2, and p is the vector of corresponding nodal forces. 

The stiffness matrix can be written as the sum of the stiffnesses due to 
normal and shear components of the elastic properties 

k = k n + k s (3.1-5) 
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3.1 

where 


(Conti nued) 





-^32*32 

x 32 2 


Symmetric 

Et 

-y32y3i 

vx 32 y 3 i 

2 

*31 

-vy 31 x 31 

-y3iy2i 

v *31 x 21 


4A(l-v 2 ) 

v -^32 x 31 

^32^21 

-vy 3 2 x 2i 

" x 32 x 31 

-vx 32 y 21 

x 32 x 21 

x 31 2 

vx 31 y 2l 

“ x 31 x 2l 


■vy 2 i x 2l X 21 ' 


and 

L_ 


_J 


x 32 2 

_x 32 y 32 

y 32 2 

Symmetri c 

y - Et 

" x 32 x 31 

y32 x 31 X 31 


s 8A(l+v) 

><32*31 

-y 32 y3i -x 3 iy 31 

y3i 2 


x 32 x 21 

-y 3 2 x 21 -X 31 X21 

y 3 l x 2! X21 2 


_2 X 32*21 

y32y 2 i x 3iy2i 

■* 31*21 - x 2iy2i *2i 2 _ 

Plane -Strain Stiffness Matrix 

- The plane-stress 

and plane-strain matrices 


of elastic coefficients are defined by Equations 2.1-4 and 2.1-5, re- 
spectively. If in the plane-stress elasticity matrix we substitute E for 
E, and Z for v, where 


E = E/(l-v 2 ) 

v = v/(l-v) 


(3.1-6a) 
(3. 1 -6b) 


we obtain the elasticity matrix for plane strain. This means that the 
element stiffness matrix for plane-strain conditions will be defined by 
Equation 3.1-5, providing that the substitutions E for E, and v for v 
are made in the matrices k n and k $ . BOPACE therefore uses the same 
subroutine to generate stiffness matrices for both plane-stress and plane 
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3.1 (Continued) 

strain problems, by using the appropriate values for Young's modulus and 
Poisson's ratio in the plane-stress matrix generation. 


Elasto-Plastic Stiffness Matrices - The complete incremental force- 
displacement relation for an element corresponds to the stress-strain 
relation 2.5-6. Because strains are constant within an element and are 
directly relatable to nodal displacements, the relations are identical 
in form: 


Ap. 


(Akf 


1J 




1J 




(3.1-7) 


where each of the force-displacement matrices k is obtained directly from 
the corresponding stress-strain matrix D, by using Equation 3.1-3. 

For example, 



A D 


Pi 


mn mi mj 


(3.1-8) 


3.2 COORDINATE SYSTEMS 

Coordinate systems are used in BOPACE to define the locations of nodes and 
to define the directions of nodal displacements. 


Location Coordinate Systems - BOPACE allows the user to define nodal 
locations by either of two global coordinate systems. These are the 
X-Y Cartesian system and the R-e cylindrical system, shown in Figure 
3.2-1. The X-Y and R-9 systems have a common origin, 0, and the angle e 
is measured counter-clockwise from the X axis. The global X-Y system 
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3.2 (Continued) 

also serves as the local element x-y system which was discussed in 
Section 3.1. 

Direction Coordinate Systems - BOPACE allows the user to define directions 
of the nodal displacements by three types of coordinate systems. These 
are local Cartesian systems, with origins at each particular node. The 
first two types are defined by the X-Y and R-e systems, respectively, 
and the third type consists of special user-defined systems. The 
three direction-coordinate types are shown in Figure 3.2-2 for a parti- 
cular node, and are denoted by the respective subscripts 1, 2 and 3. The 
first system measures displacements parallel to the global X-Y axes. 

The second system measures displacements parallel and perpendicular to 
the radius passing through the particular node. The user-defined systems 
are at given orientations, defined by the counter-cl ocRwise angle from 
the global X axis to the X-j axis. The special user-defined systems are 
useful in cases where nodal displacements are to be prescribed in 
particular directions. 




FIGURE 3.2-1: LOCATION COORDINATE FIGURE 3.2-2: DIRECTION COORDINATE 

SY5TEMS SYSTEMS 


3-6 



D5-1 7266-1 


3.3 COORDINATE-RELATED TRANSFORMATIONS 

This section gives the necessary transformations which convert BOPACE 
quantities defined in one coordinate system to the equivalent quantities 
in another coordinate system. The quantities considered are coordinates, 
displacements, forces, stiffnesses, strains and stresses. 



Basic Transformation Matrix - Referring to Figure 3.3-1, the transfor- 
mation relations will be developed between the x-y and X-Y coordinate 
systems. The angle e is the angle, measured counter-clockwise from the 
X axis to the x axis. Using the shorthand notation c and s to denote 
respectively the Cosine and Sine of the angle e (i.e., the direction 
Cosines of x given in the X-Y system), we define the 2-dimensional 
transformation matrix, r, by: 

F = [-s c] (3.3-la) 
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3.3 (Continued) 

Because of the orthogonality properties of this matrix, the inverse of 
r is equal to its transpose: 



In the following transformations, upper and lower case letters will be 
used to denote quantities which are referred to the X-Y and x-y coor- 
dinate systems, respectively. 

Coordinates - The coordinate transformations are given by 


x i 

' r U x o 


(3.3-2a) 

and 




X i 

= rT 1 . x . = 

rT. x. 

(3.3-2b) 

ij J 

ij J 

where 




X 

= {)}and X 

X >- 
!l 


Displacements - The 

displacement relations are 


*1 

= r . . q. 

ij 


(3.3-3a) 

and 




Qi 

11 

—1 

1 

Cj. — > 
Xi 

—J* 

il 

rT. q . 

(3. 3 -3b) 


where 


q - {“} and Q = {"} 
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3.3 (Continued) 

Forces - The generalized nodal force relations are 


p i = 

r ij p j 

(3.3-4a) 

and 



Pi = 

(rT >"' P j ' r ij P j 

(3.3-4b) 

Stiffnesses 

- The transformation relations for a nodal 

(2x2) stiffness 

matrix partition are 


K. . = 

k r . r . 

(3.3-5a ) 

TJ 

mn mi mj 

and 



k ij = 

'Snn r mi r mj “ 'Sun r mi r mj 

(3.3- 5b ) 


Strai ns - Transformation of the strain tensor is described by 



(3.3-6a) 
( 3 . 3- 6b ) 


where here t and E denote the strain tensor expressed in the x-y-z and 
X-V-Z coordinate systems, respectively, and r is the transformation matrix 
(of x-y-z direction cosines) which gives x-y-z coordinates in terms of 
X-Y-Z coordinates. 


The 2-dimensional (x-y plane) strain transformations are written as 



(3.3-6c) 
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3.3 (Continued) 


and 


XX j 

cc 

ss 

2cs" 

m 

X 

X 

yy £ 

ss 

cc 

-2cs 

m 

-< 

-< 

xy •• 

“CS 

cs 

cc-ss_ 

m 

X 

-< 


(3.3-6d) 


where, for example, cs denotes the product of the Cosine and Sine of 
the angle e. 


Stresses - The stresses transform between the Cartesian coordinate 
systems x-y-z and X-Y-Z in exactly the same manner as do the strains. 
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4.0 SOLUTION METHOD 

4.1 BASIC SOLUTION REQUIREMENTS 

The exact elasto-plastic-creep analysis of a structure requires the 
satisfaction, at all points in the structure, of three requirements: 

1. satisfaction of constitutive theory, which is summarized by the 
appropriate incremental stress-strain relation 

2. compatibility of strains 

3. equilibrium of stresses 

The following paragraphs summarize the BOPACE solution approach as it 
relates to satisfying these three requirements. 

Stress-Strain Relation - The incremental stress-strain relation is de- 
fined by Equation 2.5-6, and it is satisfied exactly in the BOPACE solu- 
tion procedure, provided the increment is sufficiently small so that 
incremental quantities can be treated in a differential manner. Be- 
cause of the state of constant strain and stress which exists within 
a CST finite element, the stress-strain relation for each element holds 
over the entire area of that element. 

Compatibility - Compatibility is satisfied exactly within each element 
as a result of the finite-element derivation. In the global sense, i.e., 
over the entire structure, compatibility is also satisfied exactly, 
by merging the element degrees of freedom into global degrees of freedom 
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4.1 (Continued) 

and thereby establishing the equality of displacements at appropriate 
adjacent nodes. 

Equi 1 i brium - Equilibrium is satisfied exactly within each CST element 
because of its constant stress state (as long as there are no distributed 
loads). Stresses are not necessarily in equilibrium between adjacent 
elements, although stress equilibrium is satisfied there in the limit as 
the finite-element mesh is refined. For any mesh representation of the 
structure, global equilibrium is satisfied in BOPACE in an average sense, 
because equilibrium is established between the generalized nodal forces 
defined according to the usual finite-element procedure. 

4.2 COMPARISON' OF COMMON SOLUTION METHODS 

The common stiffness methods used for solution of elasto-plastic problems 
can be classified by three general types: 

1. The pure "tangent stiffness" method 

2. The "constant-stiffness residual-load" method 

3. "Combined" methods 

Tangent-Stiffness Method - The pure tangent-stiffness method obtains 
the solution for each load increment by a single solution of the incre- 
mental equilibrium equation: 

4P i ■ 4, JOj (4.2-1) 
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4.2 (Continued) 

in which aP and aQ are the global incremental forces and displacements, 
respectively, and K 1 ^ is the Jacobian (tangent-stiffness) matrix. This 
is the type of solution used in NASTRAN's "piecewise linear analysis", 
for example. There is no equilibrium check, and no iteration is per- 
formed to improve the incremental solution. The matrix K J is determined 
by evaluation or extrapolation at previous solution points. Because in 
an actual structure the stress-strain slopes, creep rates, direction of 
the incremental plastic and creep strain vectors, etc., will generally 
vary within an increment, the pure tangent-stiffness approach can result 
in a substantial departure from the true force-displacement path unless 
load increments are kept quite small. 

Constant-Stiffness Residual-Load Method - This solution method [9] 
employs an iterative procedure. In each iteration the residual (un- 
balanced) forces are computed based on the current estimate for the 
incremental configuration, and are then applied to the constant elastic 
stiffness matrix in order to solve for displacement corrections. The 
approach is computationally efficient because it requires the formation 
and decomposition of only a single stiffness matrix, but it is not 
directly applicable to highly nonlinear structures because of con- 
vergence difficulties. 

Combined Methods - Various combined methods have been employed for 
solution of elasto-plastic problems, for example that described in 
Reference 12- These involve the use of an equilibrium check through 
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4.2 (Continued) 

the calculation of unbalanced forces, as well as various procedures for 
updating the approximate Jacobian matrix. 

BOPACE Approach - BOPACE uses a combined approach for solution, with 
the iterative procedure consisting of two basic parts: 

1. Calculation of residual nodal forces based on the estimated 
configuration and "exact" constitutive theory. 

2. Improvement of the solution configuration by reducing the residual 
forces. 

Several user controlled options are available in BOPACE for updating the 
Jacobian matrix. 

4.3 CALCULATION OF UNBALANCED FORCES 

It is assumed for the present discussion that the exact solution configura- 
tion is known at the start of a particular load increment. (Actually 
the BOPACE program takes any unbalanced forces which might remain from 
the previous increment and adds them to the present load increment, in 
order to achieve greater accuracy.) For a given iteration within the 
present increment, i.e., for a given estimate of the incremental solution, 
it is necessary to compute the corresponding unbalanced forces. This 
section describes the steps involved in computing these forces, including 
determination of element strains, element stresses, and forces. A flow- 
chart for these calculations is given in Figure 4.3-1. 
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4.3 (Continued) 

Element Strains - For the given estimate of incremental global displace' 
ments, AQ, the corresponding element displacements, Aq, are obtained by 
coordinate transformations of the type 3.3-3a. Element strains, Ac, 
are then computed using the strain-displacement matrix of Equation 3.1-2: 

^■j = r -j j Aq. (4.3-1) 

The strain Ae is the total (physical) strain increment in the element: 

Ae. = Ae® + Ae? + AeJ + Ae® (4.3-2) 

The thermal strains, Ae*, and creep strains, Ae C , are determined as des- 
cribed in Sections 2.2 and 2.4, respectively. Subtracting these strains 
from the total strain, gives: 

Ae? +P = Ae? + Ae? = Ae,. - Aef - Ae® (4.3-3) 

Element Stresses - Incremental stresses are then computed using the 
elasto-plastic relation 2.5-6 (see also Reference 13): 

A« f - AD, je f + Aef» (4.3-4) 

eO 

where e are the known cumulative elastic strains at the beginning of 
the increment. This relation holds only if the element is plastic 
throughout the load increment. Therefore, the calculation of A a may 
need to be modified, depending on which of three conditions exists 
within the particular element: 
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4.3 (Continued) 

Condition I element is elastic at end of load increment, 
i.e., either the element remains elastic or 
unloading occurs. Compute stress and elastic 
strains. Plastic strains are zero. 

Condition II element is plastic throughout load increment. 

Compute stresses, elastic and plastic strains 
from Equation 4.2-4. 

Condition III element is initially elastic, but becomes 

plastic at some point during the load incre- 
ment. Find intermediate point at which 
yielding occurs (this requires solving a simple 
quadratic equation). Compute stresses and 
elastic strains up to that point. Compute 
stresses and strains beyond yielding as for 
Condition II. 

The condition at the beginning of the increment is known for each element. 
The condition at the end of the increment is assumed, for the first 
iteration, to be Condition I. The end condition is re-evaluated during 
each iteration, using either the yield point or the plastic-strain 
vector. For an elastic element, it is determined whether or not the 
current yield point has been exceeded. For a plastic element, the 
plastic strain vector (normal to the yield surface) is observed; an out- 
ward vector (x>0) implies a plastic condition, while an inward vector 
(AiO) implies elastic unloading. 
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4.3 (Continued) 

Element Forces - The element force-stress relation for the CST element 
is directly apparent from the expressions for the strain-displacement 
and stiffness matrices given in Equations 3.1-2 and 3.1-3. The relation 
is 

Ap i = A r ij Aa j (4.3-5) 

where A is the element area, r is the strain-displacement matrix, and p^ 
are the element nodal forces. 


Global Unbalanced Forces - Incremental global forces, aP, are obtained 
from the incremental element forces, by adding nodal contributions from 
all elements and applying transformations of the type 3.3-4a. The global 
unbalanced forces, 6P, are then determined by subtracting these computed 
(internal) incremental forces from the applied (external) incremental 
loads : 

6P i = (a Load).. - AP i (4.3-6) 


4.4 IMPROVING THE SOLUTION 

The basic global relation for incremental forces and displacements cor- 
responds to the element relation 3.1-7: 


AP 


i * Kj + 4 K ij> of * - 


*?3> < p 


(4.4-1) 


where the first term accounts for forces due to change in elastic proper- 
ties, and the second term accounts for forces due to change in elastic 
strains. 
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4.4 (Continued) 

In order to improve a given displacement configuration, the displacement 
corrections 6Q corresponding to unbalanced forces 6P, are obtained in 
BOPACE by solving a set of linear equations of the form 

5P i = K-j «Qj (4.4-2) 

The matrix is the Jacobian (tangent-stiffness) matrix, or some approxi- 
mation to the Jacobian. The purpose of this section is to discuss the 
procedure for relating Equations 4.4-1 and 4.1-2, and describe BOPACE 
options for updating the Jacobian. 

Procedure - In the iterative BOPACE approach, the only global solution 
employed is the displacement-correction relation 4.4-2. The best 
approximation for the Jacobian is 

K J = K el + K pl (4.4-3) 

where K e ^ and K p ^ are evaluated at the end of the current load increment. 
The effects of change in elastic properties (/\K e and AK P ), as well as 
the effects of thermal and creep strains, are computed on the elemental 
level and accounted for by the unbalanced forces. Thus Equation 4.4-1 
is satisfied in an iterative fashion. 

Updating the Jacobian - In order to account for possible large-scale 
elastic unloading of the structure under cyclic load conditions, one or 
more initial iterations are performed for each load increment using 
only the elastic portion, K el , of the K J matrix. Succeeding iterations 
use the total matrix. 


4-9 



D5-1 7266-1 


4.4 (Continued) 

Initially the K J matrix is taken to be the usual elastic stiffness matrix 
for the structure, with elastic properties evaluated at the fabrication 
temperature. Whenever convergence is not achieved within a specified 
number of iterations, the Jacobian matrix is updated. BOPACE allows 
five options for updating the matrix K J and/or its component matrices 
K e1 and K pl . 

1. Use only elastic matrix K e ^ with no updating. This option corresponds 
to the constant-stiffness residual-load method, and is most effective 
for problems with small plastic strains and elastic properties which 
do not vary much with temperature. 

0 1 

2. Update only K . This option is best for problems with small plastic 
strains and elastic properties which vary considerably with tempera- 
ture. 

3. Update only K p ^. This option is best for problems with large plastic 
strains and elastic properties which do not vary much with temperature. 

4. Update total K^ 1 matrix, but not its component matrices. This option 
may be used for problems with large plastic strains and elastic 
properties which vary somewhat with temperature. 

5. Update both K^, and K e ^ and K p ^ matrices. This is the most effective 
option for problems with large plastic strains and elastic properties 
which vary considerably with temperature. 
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4.5 SUMMARY OF BOPACE SOLUTION METHOD 

An outline of the BOPACE solution method is given in the flowchart of 
Figure 4.5-1. In step 1, the Jacobian is initialized to the elastic 
stiffness matrix, based on elastic properties at the fabrication tem- 
perature. 

At the start of each load increment (step 2) the residual forces 6P are 
set equal to the increment of applied loads. Also, if any residual 
forces remain from the previous load increment, these are added to SP. 

The estimate for incremental displacements, AO, is set to zero. 

The iteration loop involves successive improvement of the solution, by 
solving for displacement corrections using the unbalanced forces 
and the Jacobian, and then recomputing the unbalanced forces cor- 
responding to the new displacement configuration. The displacement 
corrections 6Q are determined in step 3, and in step 4 the improved 
incremental configuration aO is updated by addition of 60. Although 
convergence of this iterative process is usually quite good, BOPACE 
has a feature for modifying the process if convergence is not occurring. 
This involves using only a specified fraction of the computed correction, 
e.g., aQ +- aQ + 0.5 6Q. This would increase the numerical stability 
but could tend to slow down convergence. 

In step 5 the strain-displacement relations are used to compute the total 
incremental strains Ae from displacements AQ. The thermal and creep 
strains are then determined in step 6, and subtracted from total strains 
to give the elastic+plastic strains required for the calculation of 
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4.5 (Continued) 

incremental stresses. Step 7 involves the major iteration algorithm, 
in which the strain is separated into elastic and plastic components. 
Incremental stresses are determined according to Equation 2.5-5, and 
the corresponding unbalanced forces are computed. 

If the maximum allowable iterations have been exceeded, step 8 is used 
to update the Jacobian matrices according to the specified updating 
option. The Jacobian update is based on the current estimates of 
the yield surface and flow parameters for each element at the end of 
the present increment. Iteration is stopped when an error norm 
(determined by a ratio of residual forces to applied forces) is 
sufficiently small. 
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5.0 LINEAR EQUATIONS 

5.1 LINEAR EQUATION FORMATION 

Single Element - Figure 5.1-1 shows a linearly elastic spring. It is 
defined by its nodal points 1 and 2, and its stiffness coefficient, k 
(force/displacement). Forces F x and displacements u are assigned to the 
nodes as shown in the figure. It is then desired to determine the stiff- 
ness matrix for this element. To do this consider all possible displace- 
ment modes which the spring is capable of experiencing. 


F X] * u ] 


O-^ 

1 


mMmm 




F X2» U 2 


Figure 5.1-1: ELASTIC SPRING AS A FINITE ELEMENT 


For example, 

a. u-j = 0, u 2 f 0 

For this case it is clear that Fx£ must be an applied load, while Fx-j is a 
reaction. Also, U 2 is an unknown displacement. From the definition of k, 

k = or Fx? = ku? 

1*2 ^ 

Equilibrium of forces gives, 

Fx-| = -Fx 2 

b. u 1 t 0, u 2 = 0 
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5.1 (Continued) 

This is the second (and final) displacement state which the spring element 
can undergo. Conditions are now reversed from those which applied in the 

first case; that is, Fx-| is now an applied load and u-j are unknown displacement. 

The applicable equations are therefore, 

k = or Fx-j = ku-| 

U 1 


and 

Fx 2 = -Fx-| 

Now assume the superposition of the individual displacement modes. Then 
/ 0 and u 2 t 0. From cases (a) and (b) above it then follows that, 


Fx-| = ku-j - ku£ 

Fx 2 = -ku-| + ku 2 
or in matrix form, 

( u d 
h) 

or simply, 

{F} = [K] {U} 



(5.1-la) 


(5.1-lb) 


In this equation [K] represents the element stiffness matrix. It is seen 
to relate nodal displacements and corresponding nodal forces. The most 
important characteristics of a finite element are represented by its 
stiffness matrix. 
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5.1 (Continued) 

Assemblage of Spring Elements - Consider now an assemblage of spring 
elements as shown in Figure 5.1-2. Nodes are given as 1, 2, 3, and 4. 

The elements are connected to each other at the nodes. Each spring has 
a different stiffness as designated by the constants k a , k^, and k c . 

It is then desired to find the stiffness equation for this assemblage. 

k fl k b k c 

12 3 4 

Figure 5.1-2: ASSEMBLAGE OF FINITE SPRING ELEMENTS 

Nodal forces are Fx^ , i = 1, 2, 3, 4 and the corresponding displacements are 
simply u.j . All possible displacement configurations for this system can be 
obtained by superimposing the four basic states. These occur as each node is 
displaced in turn, the other nodes in the meantime being held fixed. As in 
the case of the single element these displacement states can be imposed on the 
assemblage. The relevant force-displacement expressions can then be written 
for each case as follows: 

a. U-| t 0, u.j = 0, i = 2, 3, 4 

Only spring 1-2 is strained. Fx-j can be considered an applied load and Fx^ a 
reaction. No forces are carried through to nodes 3 and 4. Hence, 

k a = or Fx-j = k aU] 
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5.1 (Continued) 

and by equilibrium, 

Fx 2 = -Fx, ■ -k,u, 

b. li 2 / 0, u. = 0, i = 1, 3, 4 

Both springs 1-2 and 2-3 are now strained. Fx 2 is the applied load. Reactions 
occur at nodes 1 and 3. No load carries through to node 4. Let, 


Fx 2 - (Fx 2 ) a + ( Fx 2 ) ^ 


so that, 


k (Fx^ 

a 


then, 


Fx 2 = (k a + k b ) u 2 


k = ( px 2 )b 

K b 

u. 


while, 

Fx, ■ -(Fx 2 ), - -k a u 2 

Fx 3 = "t Fx 2^b = _l< b u 2 
c. U 3 t 0, u n - = 0, i = 1 , 2, 4 

This is similar to case (b). By a similar calculation it follows that, 

Fx 3 = (k b + k c ) u 3 

Fx 2 = -kb u 3 anc * Kx 4 = ~ k c u 3 
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5.1 (Continued) 

This is similar to case (a). Hence, 

Fx 4 = k c u 4 and Fx 3 = - k c u 4 

e. t 0, i = 1 , 2, 3, 4 

This is the general case. It is obtained by superimposing the first four 
cases. Collecting terms from the previous force-displacement equations yields, 

Fx l = k a u l - k a u 2 

Fx 2 = " k a u l + ^ k a + K b^ u 2 " k b u 3 

Fx 3 = ” k b u 2 + {k b + k c } U 3 ■ k c U 4 

Fx 4 = -k c u 3 + k c U 4 
In matrix form, 

( Fx l ) 

I Fx 2 f 

] Fx 3 ( 

( F><4 ) 



or again, 

{F> = [K] {u> 


(5.1 -2b ) 


5-5 



D5-1 7266-1 


5.1.1 (Continued) 

The assemblage stiffness matrix is shown in Equation 5.1-2a. Note that 
its order (4 x 4) is governed by the number of distinct nodal displacements 
which the assemblage can experience. Some other important features about 
[K] are as follows: (a) it is a square, symmetric matrix; (b) the stiffness 

matrix for an assemblage is only sparsely populated with non-zero elements 
and these can be banded along the main diagonal; (c) the sum of elements in 
any column (or row) is zero; (d) the stiffness matrix is singular or the 
determinant of [K] vanishes. These properties are important and except 
for (c) are generally true for any problem*. They will be referred to from 
time to time in this document. 


A useful form for writing a stiffness matrix is shown in Equation 5.1-3. 
The matrix is the same as previously given in Equation 5.1-2a. 

U 1 u 2 u 3 u 4 


[K] 


k a " k a 


“ k a ^ k a + k b^ 
0 -k k 


0 


0 

” k b 

( k b + 


-k„ 


(5.1-3) 


Note that the column of displacements { u > in Equation 5.1-2a has simply 
been turned into a horizontal position and written above the columns making 
up [K]. Actually Equation 5.1-3 is sufficient for defining the total 
stiffness equation since the order of terms in the force column {FI must 
agree with that given for {ul. See Equation 5.1-2a. 


* Vector sum of "forces" are zero; however moment terms require products 
of moment arms and forces in the moment equilibrium equations. 
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5.1 (Continued) 

Fortunately a simple, direct procedure exists for forming the assemblage 
stiffness matrix. In other words it is not necessary to use the procedure 
just described for the spring assemblage. Instead the total assemblage 
stiffness matrix can be formed directly from the stiffness matrices for the 
individual elements. 

Assemblage Stiffness by Superimposing Element Stiffnesses - It can be 
shown rigorously and is here demonstrated for the assemblage of Figure 
5 . 1-2 that superimposing the stiffnesses of the component elements 
will yield the stiffness matrix for the assemblage. It is this pro- 
cess which characterizes the direct stiffness method. 

Prior to superposition a simple preliminary step is necessary; namely, to 
increase the order of each element stiffness matrix to that applying to 
[K] for the total assemblage. This is done by simply adding columns 
of zeros (plus corresponding rows) for each displacement inapplicable to 
the given element. For example, for spring 1-2 of Figure 5.1-2, 
displacements U 3 and U 4 are irrelevant; hence, 
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5.1 (Continued) 

The columns and rows of zeros are seen to apply to and u^. Likewise 
for the other two springs. 


[K b ] 


A_ 

1 

1 

1 

1 

U 4 

0 

° ! 

0 

0 

0 

1 

1 ^ 

•cT 

1 

0 

0 

" k b 1 
1 

k b 

0 

0 

0 1 

0 

0 



u z_ 

1 

1 

^3 _ 

. A 

0 

0 

0 

0 

0 

0 

i 

0 

0 



-J- 

— — 

— 

0 

0 

1 

1 

k c 

“ k c 

0 

0 

1 

1 

" k c 

k c 


[K c ] - 


Summing the element stiffness matrices is now seen to lead directly to [K] 
for the assemblage as given by Equation 5.1-2a or Equation 5.1-3. In 
carrying out this superposition the 2x2 non-vanishing parts of the 
element stiffness matrices are seen to lie along the main diagonal as 
illustrated by Equation 5.1-4. As a result the tendency of the gross 
stiffness matrix to be sparsely populated by non-zero elements, which are 
banded along the main diagonal, can be appreciated. This fact is important 
in the numerical calculations which lead to the solution for the unknown 
displacements, particularly for large problems. 


[K] 


U] U2 U3 U4 



(5.1-4) 


5-8 



D5-1 7266-1 


5,1 (Continued) 

Shaded areas in Equation 5.1-4 represent regions where element stiffnesses 
overlap and hence jointly contribute to the gross stiffness. Non-zero 
terms lie inside the three boxes labeled Ka, Kb, and Kc respectively. 

Forming the assemblage stiffness matrix may be regarded as the key to the 
stiffness method of solution. The procedure as illustrated above for 
carrying out this step is well suited to machine operation. The 
computer determines the individual element stiffness matrices and then 
combines these to form the gross stiffness matrix. 

Special Equation Generation - Figure 5.1-3 shows two triangular plates 
connected at a single node. 



Figure 5.1-3: TWO PLATE ASSEMBLAGE 
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5.1 (Continued) 

The two plates are defined by their nodal points 1, 2, and 3, and 3, 4, and 
5. Forces Fx and Fy, and displacements u and v, are assigned to each 
node as shown for node 5 in the figure. Each plate has a 6 x 6 stiffness 
matrix as shown in Equation 5.1-5. 
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5.1 (Continued) 

Merging these matrices using superposition as described in the preceeding 
section yields a 10 x 10 matrix for the entire assemblage as shown in 
Equation 5.1-6. Note that the displacements u 3 and v 3 are constrained 
to be equal for the two elements. 

Structural discontinuity is idealized by allowing one or more of the 
displacements at a connecting node to be unconstrained in the displacement 
of the attached elements, while the remaining displacements of the node 
are constrained to be equal. 

The equation generation procedure is modified to handle this condition. 

In the problem of Figure 5.1-3 assume that the v displacements at node 3 
of the two elements are constrained to be equal, and the u displacements 
are constrained. Merging the elemental stiffness matrices of Equations 
5.1-5a and 5.1-5b for this condition yields an 11 x 11 matrix for the 
entire assemblage as shown in Equation 5.1-7. 

5.2 SOLUTION OF EQUATIONS 

Introduction - The stiffness equations can be written in matrix form 
as shown in Equation 5.2-1: 

{FI = [K] {Q} (5.2-1) 

where {F} is the force vector 

[K] is the stiffness matrix 

{Q} is the displacement vector 
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5.2 (Continued) 

The solution of Equation 5.2-1 generally involves a mixed set of known and 
unknown forces and displacements. Traditionally the approach has been to 
reorder the equations collecting the known and unknown forces and 
displacements such that Equation 5.2-1 is rewritten in partitioned matrix 
form as shown in Equation 5.2-2: 



where the subscripts are 
k = known 
u = unknown 


(5.2-2) 


While this approach is efficient and yields accurate solution, the re- 
ordering necessary is a tedious detail which should be bypassed when 
possible. A somewhat different approach is the wavefront method, which 
is efficient and as accurate as any other method. The approach is to 
decompose the n x n stiffness matrix [K], then solve for the unknown 
forces { F u ) and displacements {Q u > by a forward and backward substitution. 

The solution steps are first shown here for a system of equations involving 
lonly known forces and unknown displacements. Then the modifications are 
described for the general case of mixed known and unknown forces and 
displacements . 
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5.2 (Continued) 

Solution with Known Forces - The decomposed form of Equation 5.2-1 takes 
the form 

[U] T [D] -1 [U] {Q u > = (F k ) (5.2-3) 

where [U] is an upper triangular matrix, and [D] is a diagonal matrix whose 
elements are the diagonal elements of [U]. The solution is in three steps. 


1. Decomposition 

The elements of K are 


K ij 


i-1 

z 

k=l 


ki 


til 


kj 


u.. dtI u.. 

11 11 1J 


(5.2-4) 


Since = U.^ , the elements of U are obtained successively by 
row, as 


U ij 


K ij ' 


i-1 

z 

k=l 


-1 


U 


u, 


kk ki kj 


(5.2-5) 


2. Forward Substitution 
Let Y « D" 1 U Q 
then 


i ■ ^ u k i v k ♦ u ,i v ) 


gi vi ng 


v i ■ °;I < F i - k F , u k t v k> 


i-1 


(5.2-6) 


(5.2-7) 
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5.2 (Continued) 

3. Backward Substitution 


n _i _ 1 

Y. = E D. U. b Q. + D.! U. . Q. 
1 k=i+1 11 lk k 11 111 


giving 


Qj = y. - dt] e u ik Q k 
i i n k=i+1 ik k 


(5.2-8) 


(5.2-9) 


Solution of Mixed Set - For the general case where there is a mixed 
set of known and unknown forces and displacements, the three solution 
steps must be modified. The procedure as described is for the case in 
which a single displacement (the r^ displacement) is known, with the 
corresponding r tkl force unknown. Additional known displacements with 
corresponding unknown forces would be treated in the same manner. 


The modified form of the decomposed matrix is given by Equation 5.2-3: 


[U u ] [U lr ] [U 1n ] 

T 

~ Co,,]-' 


WiP w lr ] tu ln ] 


IQ}) 



0 C U rr^ U rnJ 


[1] 


o Cu rr ] [u r „] 

< 


>*< 

(F r H 

0 0 tu nn ] 


[D* 1 ] 
L nn J 


_ 0 0 I" «1_ 


[<V. 




where Q r , and and denote known quantities. The elements of U are 
given by equation 5.2-4 except that no contribution from [Uj is distri- 
buted to the elements of [U nn ]. Detailed steps are: 
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5.2 (Continued) 

1) Decomposition 

First rows: Compute each row (1 to r-1) according to Equation 5.2-4 

and distribute contributions to later rows. 

r th row: Compute the r^ row of U according to Equation 5.2-4 but 

do not distribute contributions to later rows. 

Last rows: Compute each row (r+1 to r) according to Equation 5.2-4 and 

distribute contributions to later rows. 


2) Forward Substitution 
First rows: Compute 


{ V 


C^l ■] ] ^ {F.| } 


using Equation 5.2-7 


(5.2-11) 


and distribute to { Y r } and f Y n > . This produces the vector 


{Y} 


< Y 1> 

w 

K> 


[U, r T (V,) 

C U ln3 <V 


r 


th 


row: 


The r 


th 


row of Equation 5.2-3 gives 



the relationship 
|< F r>- C u lr ] T{V 1 >) 


(5.2-12) 


(5.2-13) 


where {Y-j} is from Equation 5.2-11 and replaces -f F r > , and 

{Y r > is set to {Q r } and contributions are distributed forward as in 
Equation 5.2-6. 
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5.2 (Continued) 

Last rows: 

Continue forward substitution by Equation 5.2-6 to obtain 
'V ■ {( F nl-[U rm ] T fQn>-[U, rl ] T (Y 1 >| 

3. Backward Substitution 
Last rows: 

From Equation 5.2-8 

<V ■ [u 22 ]-’ (d 22 ] (»„} 

rth row: 

From Equation 5.2-12 

f F r> ■ -tu lr ] n,) + [U rr ] {Q r } ♦ [U rnl ] ro„) 

First rows: 

From Equation 5.2-8 

N,>- tU,,]-’ [D n ]|lY,J - [U 1r ] CQ r } - [U lm ] (Q n )j 


(5.2-14) 


(5.2-15) 


(5.2-16) 


(5.2-17) 
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6.0 DEFINITIONS 

This section defines symbols used in the BOPACE Theoretical Manual. Symbols used 
in Section 5 (on solution of linear equations) are not included here, because 
that section uses somewhat different notation and is considered to be self 
contained. 

Variables 

a Deviatoric stress center 

b Displacement-function coefficient 

c Kinematic hardening slope 

c,s Cosine, Sine abbreviations 

k Element stiffness matrix 

p,q Local nodal forces, displacements 

r Isotropic hardening slope 

s Deviatoric (total - hydrostatic) stress 

A 

s s - a = relative deviatoric stress 

u,v Local displacements in X,Y directions 
x,y Local Cartesian coordinates 

A Elasto-plastic hardening parameter 

C Kinematic hardening matrix 

D Elasticity matrix 

E Young's modulus; global strain 

F Yield surface function 

F c Creep hardening factor 

Kinematic hardening factor 
I Identity matrix 
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6.0 (Continued) 

K Global stiffness matrix 

K J Jacobian (tangent-stiffness) matrix 

0 Origin of coordinates 

P,Q Global nodal forces, displacements 

R Isotropic hardening matrix 

T Temperature 

W Work 

R,e Cylindrical coordinates 

a Stress center of yield surface 

e Strain 

n Direction of yield surface stress correction 

p Scalar factor for yield surface stress correction 

Y Thermal coefficient of expansion 

k Cumulative plastic hardening parameter 

k C Creep hardening parameter 

u 

k Kinematic hardening parameter 

v Stress coefficients for prescribed normal stress 

v Poisson's ratio 

A Plastic fl ow parameter 

a Stress 

r Transformation matrix 

S ubscripts 

i,j,k,«.,m,n General indices 
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6.0 (Continued) 

Superscri pts 

0 Start-of-increment quantity 

1 End-of-increment quantity 

0 Known test value 

c Creep quantity 

e Elastic quantity 

p Plastic quantity 

e+p Sum of elastic and plastic quantities 

t Thermal quantity 

Special Symbols 

6 { ) Residual (corrective) quantity 

a( ) Incremental quantity 

( ) T Matrix transpose 

( )”"* .Matrix inverse 

(-) Effective quantity 

1 Summation 

^ Proportional 

^ Non proportional 
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